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We demonstrate that nonlocally coupled limit-cycle oscillators subject to spatiotemporally white 
Gaussian noise can exhibit a noise-induced transition to turbulent states. After illustrating noise- 
induced turbulent states with numerical simulations using two representative models of limit-cycle 
\ oscillators, we develop a theory that clarifies the effective dynamical instabilities leading to the 

. turbulent behavior using a hierarchy of dynamical reduction methods. We determine the parameter 

■ region where the system can exhibit noise-induced turbulent states, which is successfully confirmed 

CN| ' by extensive numerical simulations at each level of the reduction. 
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I. INTRODUCTION 



When a dynamical system is driven by external noises, its effective dynamics generally changes P, [H, 0, 0, HI- 
We usually expect that the system tends to be more random and statistically uniform, because the noise destroys 
' the spatiotemporal structures of the system. For example, Matsumoto and Tsuda [4] found that external noise 
can destabilize chaos and produces ordered behavior in the Belousov-Zhabotinsky map, and called this phenomenon 
"noise-induced order". Shibata, Chawanya, and Kaneko [H^ studied the effect of microscopic external noise on the 
collective motion of a globally coupled map in fully desynchronized states. They demonstrated that the collective 
motion is successively simplified with the increase of external noise intensity, while without the external noise a 
^ \ macroscopic variable shows high-dimensional chaos distinguishable from random motions. 

■ In this paper, we give a counter-example to this intuition. Specifically, we demonstrate that nonlocally coupled 
\ limit-cycle oscillators can undergo a noise-induced transition from uniform states to turbulent states through effective 

dynamical instabilities, with the induced turbulent fluctuations far larger than the intensity of the driving noise. 
Our starting point is a general equation describing nonlocally coupled limit-cycle oscillators subject to spatiotem- 
, porally white Gaussian noise. We first present numerical examples using two representative models of limit-cycle 

■ oscillators, the FitzHugh-Nagumo model and the Stuart-Landau model, in order to illustrate that weak external noise 
can actually cause turbulence in such systems. To theoretically investigate this noise-induced turbulent state, we 
simplify our original equation to a Langevin phase equation by means of the phase reduction method for limit-cycle 

^ I oscillators, utilizing the fact that the external noise intensity and the coupling strength between the oscillators are 
sufficiently weak. The resulting equation describes a system of nonlocally coupled noisy phase oscillators. We then 
derive an equivalent nonlinear Fokker-Planck equation from this Langevin equation by adopting the mean-field the- 
ory, which holds exactly for our nonlocally coupled oscillators. Our nonlinear Fokker-Planck equation has a constant 
^ ■ solution corresponding to all the oscillators being in a completely desynchronized state. It is linearly stable when 
\ the external noise is sufficiently strong, but, as the noise intensity is decreased, undergoes a Hopf bifurcation at a 
certain noise intensity, giving rise to limit-cycle oscillations of the phase distribution. In the vicinity of this bifurcation 
point, we can derive a complex Ginzburg-Landau equation from the nonlinear Fokker-Planck equation by applying 
the center-manifold reduction method, which governs the small-amplitude deviation of the probability density from 
the constant solution. It is well known that the spatially uniform oscillation of the complex Ginzburg-Landau equa- 
tion becomes unstable and spatiotemporal chaos develops when the Benjamin- Feir instability condition is satisfied. 
Therefore, we expect that the Langevin phase equation and the corresponding nonlinear Fokker-Planck equation also 
exhibit spatiotemporal chaos under suitable conditions. By direct numerical simulations, we will confirm that the 
amplitude turbulence typical of the complex Ginzburg-Landau equation actually arises. In addition, we also confirm 
that the phase turbulence arises near the Benjamin- Feir criticality, which is also a hallmark of the complex Ginzburg- 
Landau equation. We then examine the situation far from the Hopf bifurcation. In our system, another smaller 
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critical noise intensity is expected to exist, above which the turbulence arises from the spatially uniform oscillation 
via a long- wave phase instability. To confirm this, we derive a Kuramoto-Sivashinsky-type equation by applying the 
phase reduction method to the spatially uniform oscillating solution of the nonlinear Fokker-Planck equation. Our 
calculation shows that the phase diffusion coefficient changes its sign from positive to negative as the noise intensity 
is increased from zero, which implies the destabilization of the spatially uniform oscillating solution. By a systematic 
numerical calculation of the phase diffusion coefficient, we determine the parameter region where our system exhibits 
a noise-induced turbulent state. From the mathematical analysis of the hierarchy of reduced equations and extensive 
numerical simulations of the dynamical equation at each level of the reduction, we will conclude that the appearance 
of the turbulence can be considered as a noise-induced transition phenomenon. 

The organization of this paper is as follows. In Sec. [Til we introduce our model and present numerical results for the 
FitzHugh-Nagumo and Stuart-Landau oscillators. In Sec. IIII) we derive Langevin and Fokker-Planck equations from 
the original equation by the phase reduction method. In Sec. |IVl we derive a complex Ginzburg-Landau equation 
by a center-manifold reduction of the Fokker-Planck equation near its Hopf bifurcation point. In Sec. [Vl we reduce 
the Fokker-Planck equation to a Kuramoto-Sivashinsky-type equation near the destabilization point of the spatially 
uniform oscillation, and draw the complete phase diagram of the noise-induced turbulence. Concluding remarks will 
be given in the final section. 

II. NOISE-INDUCED TURBULENCE IN NONLOCALLY COUPLED LIMIT-CYCLE OSCILLATORS 

In this section, we introduce a general model of nonlocally coupled noisy oscillators, and numerically demonstrate 
that the model can exhibit noise- induced turbulent states using two representative models of limit-cycle oscillators. 

A. General model 

We consider a system of nonlocally coupled limit-cycle oscillators in one-dimensional space subject to spatiotemporal 
noise. The general form of the model is given by 

/CO 
dx'G{x-x')X{x',t) + V^'n{x,t). (1) 
-CO 

Here, X{x,t) represents the state of a local limit-cycle oscillator at location x and time t. The first term on the 
right-hand side describes the dynamics of each oscillator. In the absence of the coupling and the noise, it is simply 
given by X = F{X)^ which is assumed to have a single stable limit-cycle solution. The second term describes 
the nonlocal coupling among the oscillators, where K and G{x) represent respectively the coupling matrix and the 
nonlocal coupling function. Throughout this paper, we use a simple exponential function 

G{x) = ^exp{-\x\), (2) 

which is normalized to unity in the whole space domain. The last term of Eq. ([1]) represents the external noise applied 
to each oscillator, whose intensity is controlled by the parameter a. ?7(x, t) is a spat iotempor ally white Gaussian noise 
with zero mean specified by 

{r]j (x, t)) = 0, {r]j (x, t) r]j^ {x' , t')) = 2Sj^j^S {x - x') S {t - t') , (3) 

where the subscripts j and j' denote the vector components of the noise. 

Equation ([ T]) ca n naturally be derived, for example, in the following situation. Let us consider a set of equations [l2|, 

dtX (x, t)=F{X (x, t)) + kS (x, t) + v^ry (x, t) , (4) 
TdtS (x, t) = -S (x, t) + dlS (x, t)^X (x, t) , (5) 

which describe a system of mutually interacting oscillatory elements, where the coupling between the elements is 
mediated by some substance that diffuses and decays much faster than the dynamics of the individual oscillators, e.g. 
slime molds interacting through signal molecules. By considering the r ^ limit and eliminating the fast dynamics 
of S adiabatically, we obtain Eq. ([T]) with the kernel given by Eq. ([2]) from Eqs. (|4j)-([5]). 
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FIG. 1: (Color online) FitzHugh-Nagumo oscillator with parameters e = 0.5, a = 1.0, and b = 0.3. (a): Limit-cycle orbit 
and nullclines in the X-Y plane, (b): Generalized phase y) G [0,27r], which is globally defined such that cj) = uo holds 
identically. 



B. FitzHugh-Nagumo oscillators 

As the first example, we consider the case that the local oscillator in Eq. ([1]) is given by the FitzHugh-Nagumo 
(FN) oscillator. The model is explicitly given by 

dtX{x,t) = Fx{X,Y) + KxSx{x,t) + V^vx{x,t), (6) 
dtY{x,t) = Fy{X,Y) + KYSx{x,t) + ^r]Y{x,t), (7) 

where the dynamics of each oscillator is described by 

Fx{X,Y)=€-^X-X^-Y), FY{X,Y)=aX + b, (8) 

and the nonlocal coupling term Sx{x^t) is given by 



/CO 
dx'G{x-x')X{x',t). (9) 
-co 



We consider the case that the oscillators are coupled only through the variable X. We fix the parameters of the 
oscillators as e = 0.5, a = 1.0, and b = 0.3, with which the oscillators are well in the self-oscillatory regime. 

Figure HJa) shows the limit-cycle orbit of an individual FitzHugh-Nagumo oscillator in the X-Y plane together 
with its nullclines, in the absence of the coupling and the noise. Throughout this paper, we use a generalized phase 
variable to describe the oscillator (see the next section). Figure [Tfb) shows the generalized phase variable of the 
FitzHugh-Nagumo oscillator defined on the X-Y plane, which maps the state variable (X, Y) of the oscillator to a 
single real phase variable ^ G [0, 27r]. 

We performed direct numerical simulations of Eqs. (|6])-([9]), where the oscillator fields X{x,t) and Y{x^t) are 
discretized using N = 2^^ oscillators with the spatial grid size of Ax = 0.1. See Appendix [Cl for the details of the 
numerical methods. The values of the coupling strength Kx and Ky are chosen appropriately in such a way that the 
spatially uniform oscillation is stable and the system is non-turbulent in the absence of noise. 

In order to observe the average deterministic dynamics of the system, it is desirable to filter statistical fluctuations 
due to the noise. We thus introduce a space-time dependent complex order parameter with modulus R{x^ t) and phase 
0(x,t), calculated from the oscillator phase field (t){x,t) through 

/oo 
dx' G{x- x') exp (^icf) (x', t)^ . (10) 

This order parameter represents a spatial average of the complex phase factor exp(z^) of the local oscillators over the 
coupling range. 

Results for three representative cases, denoted by FN-I, FN-H, and FN-HI, are illustrated in Fig. [2] in rows. The left 
panels display the spatiotemporal evolution of the modulus R{x^ t) of the order parameter, the middle panels display 
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FIG. 2: (Color online) Numerical simulations of nonlocally coupled noisy FitzHugh-Nagumo oscillators for three representative 
cases, (a), (b), and (c): Case FN-I (Kx = 0.02, Ky = 0.08, and a = 0.0002); (d), (e), and (f): Case FN-II {Kx = 0.02, 
Ky = 0.08, and a = 0.0006); (g), (h), and (i): Case FN-III {Kx = 0.08, Ky = 0.02, and a = 0.0006). Other parameters are 
fixed as e = 0.5, a = 1.0, and b = 0.3. The left panels (a), (d), and (g) show spatiotemporal patterns of the order parameter 
modulus; the middle panels (b), (e), and (h) show instantaneous spatial profile of the local oscillator phase; the right panels 
(c), (f), and (i) show corresponding phase portraits in the X-Y plane. The number of oscillators is = 2^° and the separation 
between neighboring oscillators is Ax = 0.1. 



how the phases of the individual oscillators are distributed in space at a given time, and the right panels display the 
snapshots of the state variables of the oscillators on the X-Y plane. 

To see what happens when we increase the noise intensity from zero, we first compare the cases FN-I (top row) and 
FN-II (middle row). The coupling parameters used in FN-I and FN-II are the same, but the noise intensity used in 
FN-II is three times stronger than that used in FN-I. In the weak noise case, FN-I, the modulus of the order parameter 
(left panel) is almost uniform in space and also constant in time. The phases of the individual oscillators (middle and 
right panels) somewhat fluctuate due to the noise, but the amplitude of the fluctuation seems to be much smaller than 
that we expect for turbulent fluctuations due to a dynamical instability of the system. In the case FN-II, the noise 
is three times stronger. Now the modulus of the order parameter exhibits quite irregular spatiotemporal behavior, 
and the amplitude of the phase fluctuations is much larger than that in the case FN-I, covering the whole range from 
to 27r. Note also that the amplitude of the phase fluctuations is far larger than the applied noise intensity, which 
indicates that they are produced by a noise-induced dynamical instability of the system. 

Only with these results, it might still be suspected that what looks like turbulence in the case FN-II might actually 
not be true turbulence, and the large fluctuations of the order parameter might simply be the result of increasing 
the noise intensity by three times. To clear up this suspicion, we compare the cases FN-II (middle row) and FN-III 
(bottom row). Now the noise intensities are the same, but the coupling parameters are slightly different. As clearly 
be seen, the order parameter observed in the case FN-III is almost uniform in space and constant in time at the same 
noise level as that used in the case FN-II. Therefore, the violent order parameter fluctuations in the case FN-II should 
not simply be statistical ones due to the noise. 
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FIG. 3: (Color online) Order parameter fluctuation (i?^) — (R)^ vs. oscillator number N for the Case FN-II (open circles) and 
the Case FN-III (crosses). (•) represents the space-time average. The solid and dashed lines have slopes and -1, respectively. 
The system size is fixed as L = NAx = 102.4. 



To distinguish the two types of order parameter fluctuations more clearly, we systematically vary the total number 
N of the oscillators in each case. Since we fix the system length, we are controlling the oscillator number density, or 
equivalent ly, the number of oscillators sitting within the coupling range. If the fluctuation of the order parameter is 
simply of statistical origin coming from the finiteness of the variance should decrease as N~'^ due to the central 
limit theorem. In contrast, if the order parameter fluctuation is due to the dynamical instability of the system, the 
same quantity should remain constant even if N is varied. 

The results are shown in Fig. O where the variance of each order parameter is plotted as a function of N in double- 
logarithmic scales for the cases FN-II and FN-III. We see that there is a vast difference in the amplitude of fluctuations, 
and, as expected, its dependence on N is clearly different between the two cases; the fluctuation amplitude is almost 
constant in the case FN-II, whereas it decreases in inverse proportion to N in the case FN-III. 

Thus, we conclude that the phase fluctuations observed in the case FN-II is not merely finite-sample statistical 
fluctuations, but they are actually turbulent fluctuations generated by the effective dynamical instability of the 
system induced by the weak external noise. 



C. Stuart-Landau oscillators 



Our second example is a system of Stuart-Landau (SL) oscillators with nonlocal coupling, which is described by 
dtW{x, = (1 + iuo) W-{1^ iP) \WfW^ KSw{x, t) + ^r]w{x, t), (11) 

and 

/oo 
dx' G {x - x')W {x' ,t) . (12) 
-CO 

The Stuart-Landau oscillator is the simplest limit-cycle oscillator derived as a normal form of the supercritical Hopf 
bifurcation [7]. Each oscillator state is now described by a complex amplitude, W . Correspondingly, the coupling term 
Sw and the noise rjw are also complex variables. We take the parameters uo and as cjq = + 1-0 and K = 0.05, and 
control the remaining parameters P and a. As in the previous FitzHugh-Nagumo case, the values of the parameter P 
is chosen in such a way that the spatially uniform oscillation is stable and the system is non-turbulent in the absence 
of noise. 

We carried out the numerical analysis of this model completely in parallel with the previous case of the FitzHugh- 
Nagumo oscillators (see Appendix [C] for the numerical methods). For the Stuart-Landau oscillator, mapping from the 
complex amplitude W to the generalized phase (j) can be analytically given as 0, [1] 

(^ = argW-/?ln|W|. (13) 

Using this definition, we calculate the order parameter of the system given by Eq. (p!Q|) . 
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FIG. 4: (Color online) Numerical simulations of nonlocally coupled noisy Stuart-Landau oscillators for three representative 
cases, (a), (b), and (c): Case SL-I (f3 = 2.6 and a = 0.0001); (d), (e), and (f): Case SL-II {(3 = 2.6 and a = 0.0006); (g), 
(h), and (i): Case SL-III (/3 = 1.6 and a = 0.0006). Other parameters are given by c^o = /3 + 1.0 and i^^ = 0.05. The left 
panels (a), (d), and (g) show the spatiotemporal patterns of the order parameter modulus; the middle panels (b), (e), and (h) 
show the instantaneous spatial profile of the local oscillator phase; the right panels (c), (f), and (i) show the corresponding 
phase portraits in the X-Y plane. The number of oscillators is = 2^° and the separation between neighboring oscillators is 
Ax = 0.1. 

Figure [H summarizes the numerical results, where three representative cases SL-I, SL-II, and SL-III are compared. 
As previous, the parameter /3 is the same for SL-I and SL-II, but the noise intensity in SL-II is six times larger than 
that in SL-I. In the case SL-I, the amplitude of the phase fluctuations is quite small, and the modulus R{x, t) of the 
order parameter is almost spatially uniform and temporally constant. In the case SL-II, the amplitude of the phase 
fluctuations becomes quite large, and the modulus R{x^ t) of the order parameter exhibits complex spatiotemporal 
dynamics. Thus, the strong turbulent fluctuation arises from the spatially uniform oscillation as the noise intensity 
is increased. However, when the parameter [3 is slightly changed, SL-III, the noise with the same intensity cannot 
induce such turbulent fluctuations. 

From the same argument as the previous FitzHugh-Nagumo case, we conclude that the strong fluctuation seen in 
the case SL-II represents a genuine turbulence of the dynamical origin induced by the external noise. 

III. REDUCTION TO NONLOCALLY COUPLED NOISY PHASE OSCILLATORS 

We have observed that nonlocally coupled limit-cycle oscillators can exhibit noise-induced turbulent states. In the 
present and the subsequent sections, we will develop a theory that explains consistently the above numerical results. 
In this section, we present our first step, namely, the derivation of a Langevin phase equation from the original 
dynamical equation for the nonlocally coupled limit cycles by means of the phase reduction method. The resulting 
equation describes a system of nonlocally coupled phase oscillators, whose validity is demonstrated numerically for 



7 



the Stuart-Landau oscillators. We also derive an equivalent nonlinear Fokker-Planck equation through the mean- field 
theory, which will be the starting point for further analysis. 



A. Phase reduction 



We apply the standard phase reduction method to our system of nonlocally coupled noisy limit-cycle oscillators, 
which derives an approximate equation consisting of only the phase variable from the original dynamical equations in 
multiple variables [7] . The phase reduction is allowed when the individual local oscillators are perturbed only slightly. 
Thus, the coupling strength and the external noise intensity should be sufficiently small. As shown in the phase 
portraits in Figs. [2] and [4] (right panels), the oscillators are always in the near vicinity of the unperturbed limit-cycle 
orbits. Therefore, the parameter values used in the previous numerical analysis satisfy the above condition. 

As we already mentioned, we use a specific definition of the phase, determined on the phase space of the limit cycle 
in such a way that (j) = uj holds identically, where uj is the natural frequency of the local oscillator. This can always 
be done by an appropriate nonlinear transformation of the phase space variables of the oscillator [g", ^ . 

Details of the derivation of a phase equation from Eq. ^ are given in Appendix El The resulting phase equation 
takes the form 

/oo 
dx' G{x-x')T{(j){x,t) -(l){x',t)) +/D^(x,t), (14) 
-CXO 

where T{(j)) represents the phase coupling function between the oscillators, D the effective noise intensity that inherits 
the effect of the noise r]{x^t) in the original equation, and ^(x,t) a real scalar spatiotemporally white Gaussian noise 
satisfying 

{^{x,t))=0, {^{x,t)^{x',t'))=25{x-x')S{t-t'). (15) 

The phase coupling function r(^) can be calculated from the dynamical equations of the coupled limit cycles. The 
effective noise intensity D can also be calculated once the parameters of the original dynamical equations are given. 
Note that the above equation describes a system of nonlocally coupled noisy phase oscillators. 



B. FitzHugh-Nagumo oscillators 



We first consider the case of the FitzHugh-Nagumo oscillators. By applying the formula developed in Appendix \A\ 
the phase coupling function can be expressed as 



27r 



dX 



KxZx (A + (/.) + KyZy (A + 4>) Xo (A + </.') , 



(16) 



where Xo(0) is the X-component of the unperturbed limit cycle, and ^x(0) and ^y(^) are the phase sensitivity 
functions of the FitzHugh-Nagumo oscillator. Though the limit-cycle solution of the FitzHugh-Nagumo model can- 
not be obtained analytically, these quantities can be calculated numerically with sufficient precision using standard 
methods ^Oj. Figure [5] displays the phase coupling function r(0) calculated for the two sets of the coupling constants 
Kx and Ky used in section [iTl For both parameter conditions, the coupling functions are the in-phase type, as is 
clear from the property dr{(j)) /d(j)\(f)=o < [7]. 

The effective noise intensity D can be expressed by the original noise intensity a and the phase sensitivity functions 
Zx and Zy as 



D = — 



1 

27r 



27T 



dX a 



(17) 



C. Stuart-Landau oscillators 



In the case of the Stuart-Landau oscillators, the limit-cycle solution and the phase sensitivity function can be 
obtained analytically. The natural frequency, the phase coupling function, and the effective noise intensity are given 
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FIG. 5: (Color online) Phase coupling functions r((/)) of the FitzHugh-Nagumo oscillators. The solid and dashed curves 
correspond to the Cases FN-I and FN-II {Kx = 0.02 and Ky = 0.08) and the Case FN-III (Kx = 0.08 and Ky = 0.02), 
respectively. The in-phase coupling condition is satisfied in both cases, i.e., ((\)) I d(\)\(^^^ < 0. 



by uj = ujQ - (3, r(0 - (j)') = -i^vTT^sin(0 - 0' + Q^), and D = a{l ^ /3^), respectively. Thus, Eq. ([T4|) takes the 
following explicit form 



/oo 
dx' G{x- x) sm{(j){x, t) - (t){x' ,t) + a) + V ^ (x, t) , 
-co 

where the parameter a is given by 

a = arg (1 + i/3) . 

By changing the time scale, we can simplify the above equation as 



dt(t){x,t) = / dx G {x - x')sm{(j){x,t) -0(x',t) + a) + V^^(x,t)^ 

J — oo 

where the rescaled noise intensity is given by 



D 



K 



(18) 



(19) 



(20) 



(21) 



To see the validity of the phase reduction, we present here results of direct numerical simulations of the Langevin 
phase equation ([18]). Figure [6] displays the numerical results, which correspond to the three cases, SL-I, SL-II, and 
SL-III, treated in Section [TTl By comparing Fig. [6] with Fig. [H we can confirm that the reduced Langevin phase 
equation nicely reproduces the behavior of the original nonlocally coupled limit-cycle oscillators. 

D. Nonlinear Fokker-Planck equation 

We now transform the Langevin phase equation (p!4|) derived above to an equivalent nonlinear Fokker-Planck 
equation, which makes the following analysis far easier. To do this, we note that the Langevin phase equation may 
be viewed as describing the dynamics of a single local oscillator driven by the nonlocal mean field of infinitely many 
other oscillators sitting within the coupling range. This means that the net coupling force experienced by this local 
oscillator is a macro- variable, and therefore its statistical fluctuation can be completely negligible. For this reason, 
the nonlocal coupling term in the phase equation can safely be replaced with its statistical average. 

Let us denote by /(^, x,t) the probability density function of the phase ^ of a single oscillator at a space-time 
point (x,t). Following the above argument, we average the coupling term in Eq. ([Tlj) by the single-oscillator phase 
distribution /(^, t). The equation then takes the form of a single-oscillator Langevin equation driven by a space-time 
dependent force V(^, x,t), that is, 



dt(l) (x, t) = V X, t) + VL) ^ (x, t) , 



(22) 
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time ^ space 

FIG. 6: (Color online) Numerical simulations of nonlocally coupled noisy phase oscillators corresponding to the nonlocally 
coupled noisy Stuart-Landau oscillators for Case SL-I [(a) and (b)], Case SL-II [(c) and (d)], and Case SL-III [(e) and (f)]. (a), 
(c), and (e): Spatiotemporal patterns of the order parameter modulus, (b), (d), and (f): Instantaneous spatial profile of the 
local oscillator phase. The number of oscillators is N = 2^° and the separation between neighboring oscillators is Ax = 0.1. 



where 



/oo n27r 
dx'G{x-x') / d(l>'T{(l>-(l)')f {(!>', x',t). 
-oo ^0 



(23) 



Note that this equation now involves only one dynamical variable, 0(x,t), as a result of statistical averaging. 

The above single-oscillator Langevin equation can be easily transformed to a single-oscillator Fokker-Planck equation 
in the form 



dt 



d 



I dx' 

I J —oo 



G{x- x') 



b'T{4>-<^')f{4>' 



f{<i>,x,t) 



D 



d^fi^,x,t) 



(24) 



Since the drift velocity itself involves the distribution /(^, x, t), this Fokker-Planck equation is nonlinear. This equation 
is the starting point for further analysis. 
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IV. AMPLITUDE EQUATION NEAR THE ONSET OF COLLECTIVE OSCILLATION 

In the following sections, we will further reduce the nonlinear Fokker-Planck equation ([24|) to analyze its dynamics 
near the destabilization points. In this section, we first apply the center-manifold reduction to Eq. near the onset 
of collective oscillations. Following Ref. pjj], we derive the complex Ginzburg-Landau equation, and then present 
results of numerical simulations that confirm the theoretical conjecture on the noise-induced turbulence. 



A. Hopf bifurcation and the complex Ginzburg-Landau equation 

The nonlinear Fokker-Planck equation ([24|) has a trivial constant solution, /(^, = l/27r, corresponding to the 
completely desynchronized state of the oscillators. When the noise intensity D is sufficiently large, this constant 
solution is stable. As D is decreased, the constant solution is destabilized via a Hopf bifurcation. 

In the previous study [11], it was shown that the nonlinear Fokker-Planck equation ([24j) has a series of critical noise 
intensities, below which the system can sustain traveling waves of various wavenumbers. Among them, the spatially 
uniform oscillating solution of the phase distribution with wavenumber zero has the largest critical noise intensity. 
Thus, when we decrease D from above, the trivial constant solution gives way to a spatially uniform oscillating 
solution at a certain critical value D = Dc. 

In the vicinity of this Hopf bifurcation point, we can derive an amplitude equation describing the slow dynamics 
of the destabilized mode. We introduce a complex amplitude A{x^t) describing the deviation of f{(j)^x^t) from the 
constant solution l/27r as 

/ (^, X, t) = — + — (a (x, t) e^^^+^^-^ + A* (x, t) e-'^'^-'^A . (25) 
27r 27r V / 

As explained in Ref. fH] , we can derive the complex Ginzburg-Landau equation for this complex amplitude 

dtA (x, t) = {Dc -D)A^ ddlA -g\A\^ A, (26) 

from the nonlinear Fokker-Planck equation by the center- manifold reduction method [7] . The parameters of the above 
complex Ginzburg-Landau equation can be expressed in terms of the Fourier components Ti of the phase coupling 
function 

CO 

T{^) = ^ r, e^'*. (27) 

^ = — CO 

As derived in Ref. [ll|, they are given by 

A = argmax^, Dc = ^-^, = -A (cj + S^Fa + Fq) , (28) 



.-XF ^ AFa (F2A + r_A) . . 

2^Fa - lUTx + u 2A 

where 5R and ^ denote the real part and the imaginary part, respectively. In what follows, we will restrict ourselves 
to the case of A = 1 and > 0, i.e., the case that the first Fourier component of the phase distribution undergoes a 
supercritical Hopf bifurcation. This is actually the case with our two representative models of limit-cycle oscillators. 
When the noise intensity D is decreased below Dc^ Eq. ([26|) starts to exhibit oscillatory behavior. In this region, 
by appropriately rescaling the variables, we obtain the standard form of the complex Ginzburg-Landau equation 

UMMMMM 

dtA{x,t) = A + {l + ici)dlA-{l + ic2)\AfA, (30) 
where the real parameters ci and C2 are given by 



We use this standard form in the following discussion. 
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It is well known that the spatially uniform oscillating solution of the complex Ginzburg-Landau equation becomes 
unstable and spatiotemporal chaos develops when the Benjamin- Feir instability condition 1 + C1C2 < is satisfied. 
Furthermore, it is also known that in the near vicinity of the Benjamin- Feir line, the modulus of the complex Ginzburg- 
Landau equation tends to be uniform, and the phase component dominates the dynamics of the system. In such a 
situation, the complex Ginzburg-Landau equation can be further reduced to the Kuramoto-Sivashinsky equation 

dte{x,t) = -C2 + (1 + ciC2)dle + (c2 - c^){d,ef - ^i^i±^a^e, (32) 

where 6(x,t) is the phase of the complex amplitude and is essentially the same as the phase of the order 

parameter defined in Eq. (p!Q|) except for the sign, as we see later. When the Benjamin-Feir instability condition 
is slightly exceeded, the phase diffusion coefficient 1 + C1C2 of this equation becomes slightly negative, leading to 
turbulent behavior. In such a case, by an appropriate rescaling of the variables, we obtain the standard form of the 
Kuramoto-Sivashinsky equation [7|] 

dtSix, t) = -die + {d,ef - dte. (33) 



B. Stuart-Landau oscillators 



For the sake of simplicity, we treat the nonlocally coupled noisy Stuart-Landau oscillators hereafter. In this case, 
as seen from Eq. ([2Q|) and Eq. ([2T]) . the phase coupling function is given by the simple sine function 

r((^) = -sin((^ + a), |a|<7r/2, (34) 

which is the in-phase type coupling. The parameters of the reduced complex Ginzburg-Landau equation can be 
calculated as 

^ cos a tana , 

De = ^, ci = -=tana, ^2 = - = -^. (35) 

By choosing the parameter a appropriately, the Benjamin-Feir instability condition can be satisfied. Thus, it was 
conjectured in Refs. [11, 12] that, since the reduced complex Ginzburg-Landau equation can exhibit turbulent behavior, 
the corresponding nonlinear Fokker-Planck equation, the Langevin phase equation, and the original nonlocally coupled 
noisy limit-cycle oscillators, could also exhibit turbulent behavior under suitable conditions. 

To examine this conjecture, we conduct systematic numerical simulations near the Hopf bifurcation curve. 
Theoretical bifurcation diagrams of the complex Ginzburg-Landau equation are plotted as a function of the pa- 
rameters a and D in Fig. [71(a), and as a function of the parameters ci and C2 in Fig. [71(b). The red solid line in 
Fig. [71(a) represents the Hopf bifurcation curve, below which the complex amplitude A{x^t) starts to oscillate. In the 
numerical simulations, we fix the noise intensity at D/Dc = 0.8, which corresponds to the black line. We choose three 
representative points on this line, indicated by P, Q, and R. The green line represents the Benjamin-Feir instability 
curve 1 + C1C2 = 0, which can be expressed as a = arctan 72 from Eq. (|35|). 

Note that the parameters ci and C2 cannot change independently in the present model. From Eq. (|35]) . only pairs of 
(ci,C2) satisfying ci = — 2c2 drawn as the straight line in Fig. [T^b) can be realized. The three representative cases P, 
Q, and R correspond to three different dynamical states of the complex Ginzburg-Landau equation, namely, spatially 
uniform oscillation, phase turbulence, and amplitude turbulence, respectively [20|. On these points, we numerically 
simulate the Langevin phase equation and the corresponding nonlinear Fokker-Planck equation (and, in some cases, 
also the complex Ginzburg-Landau equation and the Kuramoto-Sivashinsky equation). Numerical methods used 
in the simulations are summarized in Appendix \C\ In drawing the figures presented hereafter, numerical results 
obtained from the different equations are appropriately rescaled to accord with the standard form of the complex 
Ginzburg-Landau equation ([3Q|) or the Kuramoto-Sivashinsky equation ([33|) . 

Figures [8] and [9] display the results at the parameter point P, where we expect stable uniform oscillation obtained 
by direct numerical simulations of the nonlinear Fokker-Planck equation ([24|) and the Langevin phase equation (p!4|) . 
Figures[8Ka) and[8l(b) plot the temporal evolution of the order parameter phase 0(x, t) defined in Eq. (p!Q|) . As expected 
from the phase diagram of the reduced complex Ginzburg-Landau equation, both equations exhibit spatially uniform 
oscillating solutions. Small non-uniformity seen in the Langevin simulation. Fig. [81(b), is due to trivial statistical 
fluctuations. Figures [9Ka) and Fig. [9Kb) display the snapshots of the instantaneous spatial distribution of the phases. 
We can confirm that the phase distributions are spatially uniform. 

Numerical results obtained at the parameter point R, where the amplitude turbulence is expected in the reduced 
complex Ginzburg-Landau equation, are shown in Figs. [TOl and [TTl Figures [TOlf a) . \Wl h) . and [TQlf c) show spatiotemporal 
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FIG. 7: (Color online) (a): Phase diagram plotted as a function of the noise intensity D and the phase shift a. The red solid 
curve represents the Hopf bifurcation line [D — Dc = cos(a)/2). For numerical simulations near the bifurcation, the noise 
intensity is chosen such that D/Dc = 0.8 indicated by the dotted curve. The Benjamin- Feir (BF) critical line (a = arctan >/2) is 
also indicated, (b): Phase diagram plotted as a function of the parameters ci = ^d/^d = tan(a) and C2 = ^g/^g = — tan(a)/2. 
Plus (Case P: a = 0.5, D/Dc = 0.8), open circle (Case Q: a = 1.0, D/Dc = 0.8), and filled circle (Case R: a = 1.2, D/Dc = 0.8) 
correspond to spatially uniform oscillation, phase turbulence, and amplitude turbulence, respectively. The Benjamin-Feir (BF) 
line (1 + C1C2 = 0) is also indicated by the green curve. 




FIG. 8: Uniform oscillation. Parameter values are a — 0.5 and D/Dc — 0.8 (Case P). Space (horizontal)-time (vertical) plot of 
order parameter phase G(x,t) for Fokker-Planck simulation (a) and Langevin simulation using N = 2^^ oscillators (b). 




space space 

FIG. 9: Uniform oscillation. Parameter values are a = 0.5 and D/Dc — 0.8 (Case P). (a): Instantaneous spatial profile of 
phase distribution function obtained from Fokker-Planck simulation, (b): Instantaneous spatial profile of local oscillator phase 
obtained from Langevin simulation using N = 2^^ oscillators. 
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(a) (b) (c) 




FIG. 10: Amplitude turbulence. Parameter values are a = 1.2 and D jDc — 0.8 (Case R). Space (horizontal) - time (vertical) 
plot of the order parameter modulus R{x,t) is shown in rescaled units with which the complex Ginzburg-Landau equation 
takes the standard form, Eq. ()30p . whose system size is L 114. Numerical data obtained form the complex Ginzburg-Landau 
equation (a), the nonUnear Fokker-Planck equation (b), and the Langevin phase equation using N — 2^^ oscillators (c) are 
compared. 




0.0 2.0 ra¥vr;m^'Vn/- -'^i - 

wave number ^P^^® space 



FIG. 11: (Color online) Amplitude turbulence. Parameter values are a = 1.2 and D/Dc = 0.8 (Case R). (a): Spatial power 
spectrum of the order parameter as a function of the wavenumber in rescaled units with which the complex Ginzburg-Landau 
equation takes the standard form, Eq. (|3Q)) . Numerical data obtained from the complex Ginzburg-Landau equation (GL), the 
nonlinear Fokker-Planck equation (FP), and the Langevin phase equation (LE) are compared, (b): Instantaneous spatial profile 
of the phase distribution function obtained from the Fokker-Planck simulation, (c): Instantaneous spatial profile of the local 
oscillator phase obtained from the Langevin simulation using N — 2^^ oscillators. 



patterns of the modulus R{x^t) of the order parameter defined in Eq. (p!Q|) , obtained by numerical simulations of 
the complex Ginzburg-Landau equation, the nonlinear Fokker-Planck equation, and the Langevin phase equation, 
respectively. As mentioned above, the order parameter defined in Eq. (p!Q|) gives the first Fourier mode of the 
phase distribution, which corresponds to the complex conjugation of the complex amplitude for the reduced complex 
Ginzburg-Landau equation as 

R (x, t) e^®(^'') = dx' G{x- x') e^^(^''0 (36) 

) 

dx'G{x-x') / d(t)' e'"*"' f{(t)',x',t) (37) 

dx' G{x- x') A* (x^ t) e-'^^' (38) 

) 

A* {x,t)e-^^''\ (39) 



f 

J — ( 



where we used the fact that the characteristic wavelength becomes sufficiently longer than the length of the the 
nonlocal coupling near the critical point in the last approximation. The three patterns are quite similar to each other. 
Figure fTTT a) shows the spatial power spectrum of the order parameter. They are also almost identical to each other. 
Figures [TlTb) and \TT\c) show the snapshots of the phase distributions obtained from the nonlinear Fokker-Planck 
equation and the Langevin phase equation, respectively. In both figures, we can observe strongly nonuniform phase 
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FIG. 12: Phase turbulence. Parameter values are a = 1.0 and D jDc — 0.8 (Case Q). Space (horizontal) - time (vertical) plot of 
the order parameter phase gradient v{x,t) = 2dxO{x,t) is shown in rescaled units with which the Kuramoto-Sivashinsky equa- 
tion takes the standard form, Eq. (|33|) . whose system size is L 76. Numerical data obtained from the Kuramoto-Sivashinsky 
equation (a), the complex Ginzburg-Landau equation (b), and the nonlinear Fokker-Planck equation (c) are compared. 
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FIG. 13: (Color online) Phase turbulence. Parameter values are a = 1.0 and D/Dc = 0.8 (Case Q). (a): Spatial power spectrum 
of the order parameter phase gradient v{x,t) = 23x0 {x,t) is plotted as a function of the wavenumber in rescaled units with 
which the Kuramoto-Sivashinsky equation takes the standard form, Eq. ()33|) . Numerical data obtained from the Kuramoto- 
Sivashinsky equation (KS), the complex Ginzburg-Landau equation (GL), and the nonlinear Fokker-Planck equation (FP) are 
compared, (b): Instantaneous spatial profile of the phase distribution function obtained from the Fokker-Planck simulation. 



distributions due to turbulent fluctuations, which again resemble with each other. These results clearly indicate that 
the turbulent fluctuations exhibited by the three equations are generated by the same dynamical instability. 

Finally, the numerical results for the parameter Q, for which phase turbulence is expected, are shown in Figs. [12] 
and [131 Figure [12] compares the results from the Kuramoto-Sivashinsky equation, the complex Ginzburg-Landau 
equation, and the nonlinear Fokker-Planck equation, where the spatiotemporal evolution of the phase gradient v{x,t) = 
2dxQ{x^ t) of the order parameter is plotted for each equation. They are remarkably similar to each other. In Fig.[T3lfa), 
spatial power spectrum of the order parameter phase gradient is plotted for each equation. They also show excellent 
agreement with each other. Figure [T3lf b) shows a snapshot of the phase distribution function obtained from the 
nonlinear Fokker-Planck equation. As expected, long-wavelength phase fluctuations are observed. 
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FIG. 14: Schematic phase diagram of the noise-induced turbulence as a function of the phase shift a and the noise intensity D. 

These numerical results confirm the existence of the noise-induced transition from the constant solution via the 
Hopf bifurcation and the Benjamin- Feir instability. 

V. PHASE EQUATION NEAR THE DESTABILIZATION POINT OF THE SPATIALLY UNIFORM 

OSCILLATION 

In the previous section, we have analyzed the nonlinear Fokker-Planck equation near the Hopf bifurcation point 
of the constant solution using the center-manifold reduction. In this section, we investigate a different parameter 
region, where the uniformly oscillating solution of the nonlinear Fokker-Planck equation loses its stability against 
phase disturbances. As a result, a transition line to the spatiotemporal chaos will be determined, which completes 
the phase diagram of the noise-induced turbulence. 



In the previous section, it was found that the constant solution of the nonlinear Fokker-Planck equation, which 
is stable for large noise intensity, loses its stability via the supercritical Hopf bifurcation as the noise intensity is 
decreased, and the order parameter exhibits spatiotemporal chaos near the Hopf bifurcation point under suitable 
conditions. However, the spatially uniform oscillating solution of all the oscillators is stable in the absence of noise, 
since we assume the in-phase coupling in the original nonlocally coupled limit-cycle oscillators. Thus, there should 
be a lower critical noise intensity, above which the uniformly oscillating solution of the oscillators loses its stability 
and gives way to the spatiotemporal chaos. The schematic phase diagram which describes this situation is illustrated 
in Fig. [HI Below we analyze the situation where the order parameter becomes turbulent as the noise intensity is 
increased from zero. We use the method of phase reduction once again to derive a phase equation describing the 
dynamics of slowly varying wave fronts. 



We here give the outline of the analysis. Details of the analytical calculations are given in Appendix [Bl First, 
we focus on the spatially uniform oscillating solution f^iO — 0o) of the nonlinear Fokker-Planck equation, where 
= (j) — Qt. Here, Q is the frequency of the collective oscillation, and the constant 0o is the initial phase that can be 
chosen arbitrary. We then allow this phase constant 6o to be slowly space-time dependent, and denote it as 6(x,t). 
Figure [T5lf a) shows a schematic picture of a spatially uniform oscillating solution, whose order parameter phase 6o 
is constant, and Fig. [T5lfb) the case that the order parameter phase Q{x^t) is slowly varying. The phase of the 
slowly modulated oscillating solution defined here is essentially the same as the phase of the order parameter defined 
previously in Eq. (p!Q|) except for the mean drift term Qt^ as long as the spatial scale of the variation is sufficiently 
long compared to the coupling length. Thus, we use the same notation 6(x,t) for both phase variables. As shown in 



A. Lower transition line to the spatiotemporal chaos 



B. Slow phase modulation of the spatially uniform oscillation 
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FIG. 15: (Color online) Schematic phase distribution function fo{0 — G) with a spatially uniform phase, G = Go (a), and with 
a spatially slowly varying phase, G = G(x,t) (b). 




FIG. 16: (Color online) Numerically obtained spatially uniform oscillating solution (a), right zero eigenfunction (b), and left 
zero eigenfunction (c). The parameter values are a — 1.2 and D/Dc — 0.5, which give v ^ —0.371652. 

Appendix [B1 we can derive an equation for the phase 0(x,t) in a closed form using the analytical procedure given in 
Ref. 0. If we truncate the series of gradient expansion retaining the first few terms, the derived phase equation has 
the form 

dtQ (x, t) = Ddle + Jl {d^Qf - Xdle + . . (40) 

where the parameters /i, and A can be calculated using the phase coupling function r(^), the spatially uniform 
oscillating solution /o(^), and its associated left- and right eigenfunctions. Among these parameters, the phase 
diffusion coefficient v is the most important quantity. When v changes its sign from positive to negative, the spatially 
uniform oscillating solution loses its stability, and spatiotemporal chaos sets in. 

C. Phase diagram for the Stuart-Landau case 

For the case of the Stuart-Landau oscillators with the phase coupling function given in Eq. (|34|) . we numerically 
calculated the phase diffusion coefficient v at various values of the phase shift a and the noise intensity and 
determined the = curve on the a-D plane. Typical functional shapes of /o(^), uo{6)^ and Uq{6) obtained from the 
numerical calculations are illustrated in Fig. [161 

In Fig. [TTTa), dependence of the phase diffusion coefficient u on a and D is illustrated, u is positive in the red 
region, and negative in the blue region. Figure [TTlfb) represents the whole phase diagram. The blue curve represents 
u = 0. The red line represents the Hopf bifurcation curve, and the green line the Benjamin- Feir criticality near 
the Hopf bifurcation. In the parameter region surrounded by the red curve and the blue curve, the phase diffusion 
coefficient u is negative, indicating that spatiotemporal chaos appears in this region. The Benjamin- Feir criticality 
line, which we obtained in the previous section by the center-manifold reduction, is also consistent with the results of 
the present phase reduction analysis. 
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FIG. 17: (Color online) Dependence of the phase diffusion coefficient v on the phase shift a and the noise intensity D. (a): v 
was numerically evaluated in the parameter region with a G [0.0, 1.4] and D /Dc G [0.1, 0.9]. V is negative in the blue region, 
and positive in the red region, (b): The red and blue curves represent the Hopf bifurcation line {D — Dc) and the noise- induced 
transition line to turbulent state — 0). The Benjamin- Feir (BF) line is also indicated by the green line. Cases SL-I, SL-II, 
and SL-III are indicated by the triangle, the open circle, and the cross, respectively. 



nonlocally coupled noisy Umit-cycle oscillators (FitzHugh-Nagumo oscillators, Stuart-Landau oscillators), X(x,t) 

I "1st phase reduction" 
nonlocally coupled noisy phase oscillators (Langevin phase equation), (f)(x,t) 

I "mean-field theory" 

nonlinear Fokker-Planck equation, x,t) — "2nd phase reduction"^ Kuramoto-Sivashinsky equation, B(x,t) 

I "center-manifold reduction" | "near the Hopf bifurcation" 

complex Ginzburg-Landau equation, A(x,t) — "2nd phase reduction"^ Kuramoto-Sivashinsky equation, B(x,t) 



TABLE L Hierarchy of equations for noise-induced turbulence. 



Let us now look back on the numerical results of the original nonlocally coupled noisy Stuart-Landau oscillators 
(Fig.[4j) and also its reduced phase model (Fig. [6]). Three sets of the parameters used in Fig. [4] and [6l SL-I, SL-II, and 
SL-III, are also plotted on the phase diagram in Fig. [TTTb). As can be seen, only the case SL-II is in the region of 
noise-induced turbulent states. The other two cases, SL-I and SL-III, are in the parameter region where the spatially 
uniform oscillation is stable. Thus, our theory consistently explains the occurrence of noise-induced turbulence. 



VI. CONCLUDING REMARKS 



We studied a system of nonlocally coupled limit-cycle oscillators subject to spatiotemporal white Gaussian noise, 
which exhibits a remarkable phenomenon called noise-induced turbulence. We considered the case that the coupling 
and the noise are sufficiently weak, and reduced the system to nonlocally coupled noisy phase oscillators by using 
the phase reduction method. We then derived an equivalent nonlinear Fokker-Planck equation from the Langevin 
phase equation, utilizing the fact that the mean-field theory exactly holds owing to the nonlocal coupling. The 
center-manifold reduction was applied to the situation where the trivial constant solution of the nonlinear Fokker- 
Planck equation becomes unstable and starts to oscillate, and we derived the complex Ginzburg-Landau equation. 
Furthermore, we also derived the Kuramoto-Sivashinsky- type equation from the nonlinear Fokker-Planck equation 
by applying the phase reduction method once again. The hierarchy of our equations derived for the noise-induced 
turbulence is summarized in Table. H Our theoretical analysis and numerical simulations thus provide strong evidence 
for the existence of noise-induced turbulence in systems of nonlocally coupled noisy limit-cycle oscillators. 

Finally, we briefly discuss how to distinguish the noise-induced effective deterministic dynamics of the system from 
the noisy patterns when the precise dynamical equation describing the system is not available. In principle, we will 
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be able to detect the bifurcations of the effective deterministic dynamics, if appropriately filtered spatial patterns 
exhibit qualitative changes as the external noise intensity is varied. For our nonlocally coupled noisy oscillators, we 
introduced the space-time dependent complex order parameter given by Eq. (p!Q|) to observe the effective deterministic 
dynamics of the system, where the nonlocal coupling function G{x) was used as the kernel function that filters the 
noisy patterns to eliminate nonessential statistical fluctuations. For experimental systems, however, there would be 
cases that the appropriate kernel function is not known in advance and should be determined empirically. The precise 
form of the kernel function is not important, so that the essential parameter is the width of the kernel function, namely, 
the coarse-graining scale of the noisy patterns. One possible method to estimate the appropriate width of the kernel 
function is to utilize a spatial correlation function C{x) of the local oscillators at a sufficiently strong noise intensity, 
where no collective oscillations arises. The correlation length of C{x) calculated in such a regime will directly reflect 
the coupling length of the system, which gives a reasonable estimate of the appropriate kernel width. We can then 
filter the noisy patterns using a suitable localized kernel function, such as the Gaussian or the exponential kernel, 
whose width is given by the obtained correlation length. For our nonlocally coupled oscillators, the correlation length 
of C{x) directly reflects the coupling length of G{x). Using the method explained above, we could actually observe 
spatiotemporal patterns similar to those observed using the complex order parameter, Eq. (p!Q|) . without assuming 
prior knowledge of the nonlocal coupling function G{x). 

We believe that noise-induced turbulence can be experimentally realized in the near future. 
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APPENDIX A: PHASE REDUCTION OF COUPLED NOISY LIMIT-CYCLE OSCILLATORS 

We consider the following Langevin equation describing an ensemble of coupled identical limit-cycle oscillators 
subject to independent noises 

dtXr^ {t) = F (X,) + ^n^n' (^n, ^nO + V^ry, (t) , (Al) 

n' 

where Xnit) represents the state of the n-th oscillator at time t, F{Xn) the individual dynamics, Vn,n'{Xn^ Xn') the 
interaction between n- and n'-th oscillators, a the noise intensity, and r]^{t) the external noise added independently 
to each oscillator. The noise ?7^(t) is assumed to be white Gaussian, whose statistics are given by 

{Vi it)) = 0, {ri {t) vC it')) = 25,j,Sn,n'5 {t - t') , (A2) 

where the superscripts j and / denote the vector component. When the coupling term and the noise term are 
sufficiently small, the phase reduction method is applicable. Applying the phase reduction method Q, we obtain the 
following Langevin equation for the phase variables: 

dt^n {t)=U^YZ • Vn,n^ ((/>n, M + V^^ (</>n) ' Vn (^) • (A3) 

n' 

Here, uj is the natural frequency of the oscillators, Z{(j)n) is the phase sensitivity function p?], and Vn,n'(0n5 0nO 
is the abbreviation of V n,n'{XQ{(j)n)^ XQ{(j)ri'))i where Xq is the unperturbed limit-cycle solution. By virtue of the 
independent Gaussian statistics, we can rewrite the above equation ()A3p in the following form 

dt^n (t) = + ^ Z {(^n) ■ Vn^n' {^n. M + 11^ (</>n)|| (t) , (A4) 
n' 

where ||Z(^^)|| = ^y Z{<pri) ' ^{(t^n) and the statistics of the white Gaussian noise is given by 

{in (t)) = 0, {in {t) {t')) = 2Sn,n'S {t - t') . (A5) 

We introduce the new slow phase variables ^n as 

(j)n=UJt^^n, (A6) 
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and rewrite the above Langevin equation (|A4p as 



(A7) 



The corresponding Fokker-Planck equation describing the evolution of the probabihty density function P{{(j)n}^t) of 
the phase variables is given by 



(A8) 



where 



^ Z (cjt + (/)n) • Vn,n' {^t ^ (t)n,(^t ^ j ■ ^ 



1 d 



((7||Z(wi + (/)„)f ) 



P{{4>n},t) 



(A9) 



Rapidly oscillating quantities in Jn({0n},^) are now time averaged over one period of the oscillator. The resulting 
averaged Fokker-Planck equation has the following form 

^P({0„},t) =^(^-^^ i^Vn,n'{^n-^n')P{{K},t)-D^j{{^n},t)^, (AlO) 



where the phase couphng function r„,„/(0„ — 4>n') is given by 



1 f'^'" 
= — dXZiX 



n' I 



and also the effective noise intensity D is given by 



D = ^ [ dXaZ{\ + ^„)-Z{\ + 4>n) . 
Jo 

The Langevin equation corresponding to the averaged Fokker-Planck equation (jAlOp is expressed as 

n' 

Finally, we obtain the following Langevin equation for the phases: 

dt(pn {t)=U) + Y^ r„,n' {<Pn - 4>n') + VD (t) ■ 



(All) 



(A12) 



(A13) 



(A14) 



Now let us consider the nonlocally coupled oscillator system. The above formulae can be written in the form 

J2 Vn,n' {Xn, X„,)^ J dx' G {x - x') V {X , X') , (A15) 

rn,n' {4>n "</-«') ^ j dx' G {x - x') F ((/. - (^') , (A16) 

where 

V{(^-(^') = ^ I Z (A + (/)) • V (A + (/), A + . (A17) 
Here, X' and (j)' are the abbreviations of X{x' ^t) and (j){x\t), respectively. 
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APPENDIX B: SECOND PHASE REDUCTION OF THE NONLINEAR FOKKER-PLANCK EQUATION 



In this appendix, we present analytical procedures to derive a phase equation describing the slow phase modulation of 
the spatially uniform oscillating solution. We assume long-wavelength modulation to the spatially uniform oscillation, 
so that the spatial derivatives of the function are small quantities. We thus expand the nonlocal coupling 

term of the Fokker-Planck equation as 



/oo 
dx' G{X- X') f (0, X', t)=J2 dl'^f {4>, X, t) , 

-oo 



(Bl) 



where G2m is the 2m-th moment of G{x). It is given by G2m = dx G{x)x'^"^ / {2my., where Go = 1 holds from the 
normalization condition. Due to the spatial reflection symmetry, only even moments G2m remain. For the coupling 
function used in this paper, G{x) = exp(— |a;|)/2, we obtain G2m = 1 for all m. The nonlinear Fokker-Planck equation 
is expanded as 



df{(j),x,t) 
dt 



{^ + J^ ^ d^'ri^-^')f{^',x,t)^fi^,x,t) 
d4>' T{4>- 4>') {G2dlf (0', X, t)}f (</., X, t) 

d(t>' r ((A - 0') {Gidtf (0', x,t)]f {cp, X, t) 



D 



d^f{<t>,x,t) 



d_ 

d(l> 
d_ 
d4> 
d 



(B2) 

Let us denote by f{(j)^x^t) = /o(^) = /o(0 — ^t) the spatially uniform oscillating solution of the nonlinear Fokker- 
Planck equation, where O is the collective frequency. Inserting this expression in the nonlinear Fokker-Planck equation, 
we find that /o(^) satisfies the following equation 



9o{9)fo{9) 



0, 



where 



go{e)= r d9'T{9-e')fo{9'). 

Jo 



(B3) 



(B4) 



Let u{9^t) represent disturbance to the spatially uniform oscillating solution defined by f{(p^x^t) = fo{0) -\-u{0^t). 
Eq. (|B2p is linearized in u{0,t), i.e., dfU = Lu, where the linearized operator L is defined as 



go {e)u{9) 



d_ 



2lT 



fo{9) d9'T{9-9')u{9') 



whose adjoint operator L*, defined by J^^ dO u* {0)Lu{0) = J^^ dO u{0)L*u* {0), is expressed as 

= D^^* {9) -{n-co) ^u* {9) + go (9) ^u* {9) + d9' T {9' - 9) fo {9') -^u* {9') . 

The left- and the right eigenfunctions and their eigenvalues of these linear operators are determined from 

Lui = Xiui, L'^ul = Xiul {I = 0, 1, 2, • • •) . 
The eigenfunctions are assumed to be orthonormalized as 

deu^iO) {9)=5im. 
Here we should note that the right zero eigenfunction can be chosen as 

Luo = 0, uo (9) = ^/o {9) , Ao = 0, 



(B5) 



(B6) 



(B7) 



(B8) 



(B9) 
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which fohows from the differentiation of Eq. ()B3p with respect to 0. 

We now apply the second-order phase reduction to the nonhnear Fokker-Planck equation (|B2p by treating spatial 
derivatives as perturbations. Following the procedure developed in Ref. 0, the Kuramoto-Sivashinsky-type phase 
equation describing the slowly varying phase modulation can be derived in the form 



where 



dtQ (x, t) = vdlQ + Jl (S^O)^ - Xd^Q 

p27T 

P = -G2 

Jo 



K iv) [ao /o iv) 



bo fo 



(BIO) 
(Bll) 
(B12) 



and 



X = G. 



Jo ' 



K {^) [ao {v) fo if ) 



d 



d^P u*i {if) — [ao {if) fo (ip) 



Here, the quantities ai{ip) and 6/(99) are defined by 



2lT 



dip'T{^-^')ui (^') 



27r 



(B13) 

(B14) 
(B15) 



In general, the spatially uniform oscillating solution fo{0) cannot be obtained analytically. Thus, fo{0) and the 
associated zero eigenfunctions should be calculated numerically in order to evaluate the phase diffusion coefficient u 
(jBlip . This can be done with sufficient precision by applying the numerical relaxation method using M = 2^ modes 
for the phase. See Appendix [Q for the details of the numerical methods. 

APPENDIX C: NUMERICAL METHODS 

We used spatially periodic boundary conditions in all the numerical simulations throughout this paper. We con- 
firmed our numerical simulation results are not changed even if we further increase the number of grid points N for 
the space or the number of modes M for the phase. 



1. Algorithm for the nonlocally coupled noisy limit-cycle oscillators 

We used an explicit Euler scheme with a time step At = 0.01 for the equation 

/CO 
dx'G{x-x')X{x',t) + ^/^r1{x,t). (CI) 
-OO 

The system is discretized using N = 2^^ grid points. The nonlocal coupling term can be efficiently calculated by using 
the fast Fourier transform (FFT) technique because it is simply a convolution form with the kernel given by Eq. ([2]). 



2. Algorithm for the nonlocally coupled noisy phase oscillators 

The system is discretized using N = 2^^ grid points for Fig. [6] and N = 2^^ grid points for Figs. ISlfTTl We used an 
explicit Euler scheme with a time step At = 0.01 for the equation in the following form 



dfcj) {x,t) = uo - R (x, t) sin(^ (x, t) - Q (x, t) + a) + ^/D ^ (x, t) . 



(C2) 
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R (x, t) exp(^ie (x, t)) = dx' G{x- x') exp(^i0 (x', t)^ , (C3) 
where the FFT is used for the calculation of the spatial convolution. 

3. Algorithm for the nonlinear Fokker-Planck equations 

We used a pseudo-spectral method with M = 2^ modes for the phase component, and spatial discretization with 
TV = 2^^ grid points for the equation 



[uj-R (x, t) sin(0 - e (x, t) + a) }/ (0, X, t)] + ^ (C4) 



where 

(x, t) exp(^ie (x, t)) = dx'G{x- x') y e^^7 (0', x', t) . (C5) 

A modified Euler scheme with a time step At = 0.01 is used for the temporal integration. FFT is used for the 
calculation of the spatial convolution. 

4. Algorithm for the complex Ginzburg-Landau equation 

We used a pseudo-spectral method for the complex amplitude using N = 2^^ modes for the equation 

dtA{x, t) = A + (1 + ici)dlA - (1 + ic2)\A\^A, (C6) 
where a modified fourth-order Runge-Kutta scheme with a time step Z\t = 0.01 is used for the temporal integration. 

5. Algorithm for the Kuramoto-Sivashinsky equation 

We used a pseudo-spectral method using TV = 2^° modes for the equation in the following form 

dtv{x, t) = -dlv ^vd:,v- dlv, (C7) 

where v{x,t) = 29^6(x,t). Temporal integration was done by a modified fourth-order Runge-Kutta scheme with a 
time step At = 0.01. 

6. Algorithm for the numerical relaxation method 

To determine the spatially uniform oscillating solution /o(^), we numerically evolve the following nonlinear Fokker- 
Planck equation 



dt 



L0 + d4>'T{<p-<p')fo{4>',t)^fo{4>,t) 



from an appropriate initial condition, which converges to a steadily rotating wave packet with the collective drift 
velocity Cl, fo {0 = (j) — Qt). The right eigenfunction uq{0) is simply obtained by diff"erentiating /o (0) by 0. To obtain 
the left eigenfunction Uq{0)^ we evolve the following equation for UQ{6,t)^ 

from an appropriate initial condition using the fo{0) and the Q obtained above. We constantly rescale u^{0^t) 
throughout the numerical evolution, so that the normalization condition J^^ dOuQ{0,t)uo{0) = 1 is always satisfied. 
After sufficient relaxation, UQ{0,t) converges to the desired left eigenfunction Uq{0). For the numerical calculations. 
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we used a pseudo-spectral method using M = 2^ modes, and a modified Euler scheme with a time step Z\t = 0.01 for 
the temporal integration. 
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